Name: Andy MacLachlan

Originality declaration

I, \[**insert your name**\], confirm that the work presented in this assessment is my own. Where information has been derived from other sources, I confirm that this has been indicated in the work.

date: 16 December, 2021

Start your response here

Initial project scope

The codes are as follow:

# load the package
library(spatstat)

## Loading required package: spatstat.data

## Loading required package: spatstat.geom

## spatstat.geom 2.3-0

## Loading required package: spatstat.core

## Loading required package: nlme

## Loading required package: rpart

## spatstat.core 2.3-1

## Loading required package: spatstat.linnet

## spatstat.linnet 2.3-0

## 
## spatstat 2.2-0       (nickname: 'That's not important right now') 
## For an introduction to spatstat, type 'beginner'

library(here)

## here() starts at /Users/mwanghoffin/Documents/GitHub/CASA0005/Practice Exam 1

library(sp)
library(rgeos)

## rgeos version: 0.5-8, (SVN revision 679)
##  GEOS runtime version: 3.9.1-CAPI-1.14.2 
##  Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
##  GEOS using OverlayNG
##  Linking to sp version: 1.4-5 
##  Polygon checking: TRUE

library(maptools)

## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.

library(GISTools)

## Loading required package: RColorBrewer

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:spatstat.geom':
## 
##     area

library(tmap)

## Registered S3 methods overwritten by 'stars':
##   method             from
##   st_bbox.SpatRaster sf  
##   st_crs.SpatRaster  sf

library(sf)

## Linking to GEOS 3.9.1, GDAL 3.2.3, PROJ 7.2.1

library(geojson)

## 
## Attaching package: 'geojson'

## The following object is masked from 'package:graphics':
## 
##     polygon

library(geojsonio)

## Registered S3 method overwritten by 'geojsonsf':
##   method        from   
##   print.geojson geojson

## 
## Attaching package: 'geojsonio'

## The following object is masked from 'package:base':
## 
##     pretty

library(tmaptools)

Set up the data

# first, get the Boundaries of Community Districts in New York City
NYDistricts <- st_read(here::here("Community Districts", "geo_export_3ed95fd1-bcd5-42fd-8e67-cf62e9e8398b.shp")) %>% 
  st_transform(., 32618)

## Reading layer `geo_export_3ed95fd1-bcd5-42fd-8e67-cf62e9e8398b' from data source `/Users/mwanghoffin/Documents/GitHub/CASA0005/Practice Exam 1/Community Districts/geo_export_3ed95fd1-bcd5-42fd-8e67-cf62e9e8398b.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 71 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.25559 ymin: 40.49613 xmax: -73.70001 ymax: 40.91553
## Geodetic CRS:  WGS84(DD)

NYDistricts

## Simple feature collection with 71 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 563069.7 ymin: 4483098 xmax: 609762.3 ymax: 4529952
## Projected CRS: WGS 84 / UTM zone 18N
## First 10 features:
##    boro_cd shape_area shape_leng                       geometry
## 1      206   42664312   35875.71 MULTIPOLYGON (((595104.1 45...
## 2      404   65739662   37018.37 MULTIPOLYGON (((597308.6 45...
## 3      203   44796871   33500.07 MULTIPOLYGON (((594364.8 45...
## 4      304   56662613   37007.81 MULTIPOLYGON (((593253.9 45...
## 5      205   38316975   29443.05 MULTIPOLYGON (((593432.1 45...
## 6      312   99525501   52245.83 MULTIPOLYGON (((586966.1 45...
## 7      483  191997989  106865.65 MULTIPOLYGON (((605960.8 44...
## 8      411  260363028  103434.76 MULTIPOLYGON (((605933.1 45...
## 9      207   53311689   44812.15 MULTIPOLYGON (((594782.6 45...
## 10     208   92071724   47817.43 MULTIPOLYGON (((592919.8 45...

qtm(NYDistricts)

# Now lets get the location of the evictions
library(tidyverse)

## Registered S3 method overwritten by 'cli':
##   method     from         
##   print.boxx spatstat.geom

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::collapse() masks nlme::collapse()
## x dplyr::filter()   masks stats::filter()
## x dplyr::lag()      masks stats::lag()
## x dplyr::select()   masks MASS::select()

library(janitor)

## 
## Attaching package: 'janitor'

## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

Evictions <- read_csv(here::here("Evictions.csv"),
                      na = c("", "NA", "n/a"), 
                      locale = locale(encoding = 'Latin1'), 
                      col_names = TRUE) %>% 
  na.omit(Evictions)

## Rows: 66559 Columns: 20

## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (12): Court Index Number, Docket Number, Eviction Address, Eviction Apar...
## dbl  (8): Eviction Postcode, Latitude, Longitude, Community Board, Council D...

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

#NewEvictions <- na.omit(Evictions)

NewEvictions <- Evictions %>% 
  st_as_sf(., coords = c("Longitude","Latitude"),
           crs = 4326) %>%
  st_transform(., 32618) %>% 
  clean_names()

summary(NewEvictions)

##  court_index_number docket_number      eviction_address  
##  Length:52917       Length:52917       Length:52917      
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character  
##                                                          
##                                                          
##                                                          
##  eviction_apartment_number executed_date      marshal_first_name
##  Length:52917              Length:52917       Length:52917      
##  Class :character          Class :character   Class :character  
##  Mode  :character          Mode  :character   Mode  :character  
##                                                                 
##                                                                 
##                                                                 
##  marshal_last_name  residential_commercial   borough          eviction_postcode
##  Length:52917       Length:52917           Length:52917       Min.   :10000    
##  Class :character   Class :character       Class :character   1st Qu.:10453    
##  Mode  :character   Mode  :character       Mode  :character   Median :10468    
##                                                               Mean   :10739    
##                                                               3rd Qu.:11221    
##                                                               Max.   :11698    
##   ejectment         eviction_legal_possession community_board  council_district
##  Length:52917       Length:52917              Min.   : 1.000   Min.   : 1.00   
##  Class :character   Class :character          1st Qu.: 4.000   1st Qu.:12.00   
##  Mode  :character   Mode  :character          Median : 7.000   Median :17.00   
##                                               Mean   : 7.899   Mean   :22.91   
##                                               3rd Qu.:12.000   3rd Qu.:36.00   
##                                               Max.   :18.000   Max.   :51.00   
##   census_tract         bin               bbl                nta           
##  Min.   :     1   Min.   :1000000   Min.   :1.000e+09   Length:52917      
##  1st Qu.:   195   1st Qu.:2008102   1st Qu.:2.027e+09   Class :character  
##  Median :   371   Median :2100990   Median :2.046e+09   Mode  :character  
##  Mean   :  8546   Mean   :2578910   Mean   :2.508e+09                     
##  3rd Qu.:   988   3rd Qu.:3216952   3rd Qu.:3.051e+09                     
##  Max.   :157901   Max.   :5171959   Max.   :5.080e+09                     
##           geometry    
##  POINT        :52917  
##  epsg:32618   :    0  
##  +proj=utm ...:    0  
##                       
##                       
## 

#plot the evictions in the city
tmap_mode("plot")

## tmap mode set to plotting

tm_shape(NYDistricts) +
  tm_polygons(col = NA, alpha = 0.5) +
tm_shape(NewEvictions) +
  tm_dots(col = "blue")

This part of R might not useful.

# join the New York District data and eviction data
#Joinfun <- function(data1, data2){

#output<- data1%>%
#  st_join(NYDistricts,.)%>%
#   add_count(boro_cd, name="eviction-borough") 

#  return(output)
#}

#NYEvictions <- Joinfun(NewEvictions, NYDistricts)

#NYEvictions

Data cleaning

This two part cannot be loaded.

#remove duplicates
#library(tidyverse)

#library(sf)
#NYEvictions <- distinct(NYEvictions)

Spatial subsetting

#NYEvictionsSub <- NYEvictions[NYDistricts,]
#check to see that they've been removed
#tmap_mode("plot")
#tm_shape(NYDistricts) +
#  tm_polygons(col = NA, alpha = 0.5) +
#tm_shape(NYEvictionsSub) +
#  tm_dots(col = "blue")

Try only extract one district

#extract the district

CD502 <- NYDistricts %>%
  filter(., boro_cd == 502)

#Check to see that the correct borough has been pulled out
tm_shape(CD502) +
  tm_polygons(col = NA, alpha = 0.5)

summary(CD502)

##     boro_cd      shape_area          shape_leng              geometry
##  Min.   :502   Min.   :591533049   Min.   :142729   MULTIPOLYGON :1  
##  1st Qu.:502   1st Qu.:591533049   1st Qu.:142729   epsg:32618   :0  
##  Median :502   Median :591533049   Median :142729   +proj=utm ...:0  
##  Mean   :502   Mean   :591533049   Mean   :142729                    
##  3rd Qu.:502   3rd Qu.:591533049   3rd Qu.:142729                    
##  Max.   :502   Max.   :591533049   Max.   :142729

Clip the eviction data to the single borough

#clip the data to our single borough
NewEvictionsCD <- NewEvictions[CD502,]

#check that it's worked
tmap_mode("plot")

## tmap mode set to plotting

tm_shape(CD502) +
  tm_polygons(col = NA, alpha = 0.5) +
tm_shape(NewEvictionsCD) +
  tm_dots(col = "blue")

The first thing we need to do is create an observation window for to carry out its analysis within — we’ll set this to the extent of the CD503 boundary

#now set a window as the borough boundary
window <- as.owin(CD502)
plot(window)

To use spatstat, we need point pattern analysis, we need to create a point pattern (ppp) object.

#create a sp object
NewEvictionsCD <- NewEvictionsCD %>%
  as(., 'Spatial')
#create a ppp object
NewEvictionsCD.ppp <- ppp(x = NewEvictionsCD@coords[,1],
                        y = NewEvictionsCD@coords[,2],
                          window=window)

## Warning: data contain duplicated points

# just to check
NewEvictionsCD@coords[,1]

##   [1] 573358.6 578156.8 568457.7 575043.7 576533.2 573475.7 571264.8 578739.8
##   [9] 573445.4 579349.6 571217.3 571217.3 579368.4 570707.0 577399.5 574357.7
##  [17] 579349.6 578940.5 577432.9 578638.3 576131.2 577828.3 578156.8 577999.9
##  [25] 570766.8 571591.3 578483.1 579325.3 573570.1 577271.1 575865.3 575512.0
##  [33] 574361.9 575330.0 568457.7 577218.0 579349.6 573191.5 571253.7 570630.6
##  [41] 570895.6 576826.6 570950.7 578959.5 570865.5 576283.9 575507.1 577543.6
##  [49] 579368.4 579205.0 578813.4 575512.0 578887.1 571622.2 579323.1 579090.7
##  [57] 570613.7 579289.4 578156.8 571130.6 571396.9 576076.3 575889.1 576106.0
##  [65] 576885.2 578461.4 570907.8 571468.0 570865.5 574357.7 578478.0 575053.3
##  [73] 575743.3 570486.8 570613.7 574577.3 571026.6 575043.7 571372.9 576077.3
##  [81] 573156.1 578156.8 576062.9 571889.8 578813.4 575878.2 574364.0 578941.7
##  [89] 574433.5 571396.9 575507.1 576471.6 578910.9 576543.9 572858.7 578481.6
##  [97] 571117.6 568559.0 575865.3 571832.7 571073.3 575865.3 576283.9 574364.0
## [105] 576077.4 568728.8 575707.5 577857.6 578057.7 574353.5 570970.1 577949.3
## [113] 575507.1 571199.3 576077.3 575512.0

# check the new ppp object
NewEvictionsCD.ppp %>%
  plot(.,pch=16,cex=0.5, 
       main="Evictions in CD502")

Data Analysis

Kernel Density Estimation

The size and shape of the Kernel affects the density pattern produced

NewEvictionsCD.ppp %>%
  density(., sigma = 500) %>%
  plot()

Ripley’s K

K <- NewEvictionsCD.ppp %>%
  Kest(., correction="border") %>%
  plot()

Kval <- as.data.frame(Kest(NewEvictionsCD.ppp, correction = "border"))

Kval

##               r         theo     border
## 1      0.000000 0.000000e+00   244878.6
## 2      4.026481 5.093323e+01   261203.8
## 3      8.052962 2.037329e+02   285691.7
## 4     12.079443 4.583991e+02   285691.7
## 5     16.105925 8.149317e+02   302016.9
## 6     20.132406 1.273331e+03   302016.9
## 7     24.158887 1.833596e+03   318342.2
## 8     28.185368 2.495728e+03   383643.1
## 9     32.211849 3.259727e+03   391805.7
## 10    36.238330 4.125592e+03   391805.7
## 11    40.264812 5.093323e+03   398679.5
## 12    44.291293 6.162921e+03   406985.4
## 13    48.317774 7.334385e+03   435725.0
## 14    52.344255 8.607716e+03   452483.6
## 15    56.370736 9.982913e+03   452483.6
## 16    60.397217 1.145998e+04   413177.0
## 17    64.423698 1.303891e+04   429186.0
## 18    68.450180 1.471970e+04   446633.9
## 19    72.476661 1.650237e+04   468923.1
## 20    76.503142 1.838690e+04   523029.6
## 21    80.529623 2.037329e+04   523029.6
## 22    84.556104 2.246155e+04   556978.8
## 23    88.582585 2.465168e+04   556978.8
## 24    92.609067 2.694368e+04   556978.8
## 25    96.635548 2.933754e+04   556978.8
## 26   100.662029 3.183327e+04   556978.8
## 27   104.688510 3.443086e+04   594110.7
## 28   108.714991 3.713032e+04   686940.5
## 29   112.741472 3.993165e+04   696223.4
## 30   116.767953 4.283485e+04   705506.4
## 31   120.794435 4.583991e+04   724072.4
## 32   124.820916 4.894683e+04   724072.4
## 33   128.847397 5.215563e+04   724072.4
## 34   132.873878 5.546629e+04   721866.5
## 35   136.900359 5.887881e+04   700284.8
## 36   140.926840 6.239321e+04   700284.8
## 37   144.953321 6.600947e+04   707656.2
## 38   148.979803 6.972759e+04   717623.2
## 39   153.006284 7.354758e+04   717623.2
## 40   157.032765 7.746944e+04   755476.5
## 41   161.059246 8.149317e+04   765549.5
## 42   165.085727 8.561876e+04   775622.5
## 43   169.112208 8.984622e+04   795768.6
## 44   173.138690 9.417554e+04   807921.9
## 45   177.165171 9.860673e+04   828505.9
## 46   181.191652 1.031398e+05   837610.4
## 47   185.218133 1.077747e+05   848015.5
## 48   189.244614 1.125115e+05   879230.8
## 49   193.271095 1.173502e+05   910041.4
## 50   197.297576 1.222907e+05   931082.8
## 51   201.324058 1.273331e+05   941603.5
## 52   205.350539 1.324773e+05   973165.7
## 53   209.377020 1.377235e+05   984100.1
## 54   213.403501 1.430714e+05  1021336.3
## 55   217.429982 1.485213e+05   979514.4
## 56   221.456463 1.540730e+05  1033931.8
## 57   225.482945 1.597266e+05  1033931.8
## 58   229.509426 1.654821e+05  1056964.3
## 59   233.535907 1.713394e+05  1101004.5
## 60   237.562388 1.772986e+05  1145044.7
## 61   241.588869 1.833596e+05  1145044.7
## 62   245.615350 1.895225e+05  1189084.9
## 63   249.641831 1.957873e+05  1244135.1
## 64   253.668313 2.021540e+05  1255145.2
## 65   257.694794 2.086225e+05  1277165.2
## 66   261.721275 2.151929e+05  1297760.5
## 67   265.747756 2.218651e+05  1297760.5
## 68   269.774237 2.286393e+05  1297760.5
## 69   273.800718 2.355153e+05  1347888.6
## 70   277.827200 2.424931e+05  1347888.6
## 71   281.853681 2.495728e+05  1347888.6
## 72   285.880162 2.567544e+05  1380367.8
## 73   289.906643 2.640379e+05  1344312.9
## 74   293.933124 2.714232e+05  1356002.6
## 75   297.959605 2.789104e+05  1356002.6
## 76   301.986086 2.864994e+05  1373537.1
## 77   306.012568 2.941903e+05  1408606.2
## 78   310.039049 3.019831e+05  1420295.8
## 79   314.065530 3.098778e+05  1443675.2
## 80   318.092011 3.178743e+05  1478744.2
## 81   322.118492 3.259727e+05  1525502.9
## 82   326.144973 3.341729e+05  1550489.6
## 83   330.171455 3.424750e+05  1550489.6
## 84   334.197936 3.508790e+05  1562325.4
## 85   338.224417 3.593849e+05  1597832.8
## 86   342.250898 3.679926e+05  1633340.2
## 87   346.277379 3.767022e+05  1602385.0
## 88   350.303860 3.855136e+05  1602385.0
## 89   354.330341 3.944269e+05  1626663.6
## 90   358.356823 4.034421e+05  1650942.2
## 91   362.383304 4.125592e+05  1687360.0
## 92   366.409785 4.217781e+05  1687360.0
## 93   370.436266 4.310989e+05  1699499.3
## 94   374.462747 4.405215e+05  1699499.3
## 95   378.489228 4.500460e+05  1735917.1
## 96   382.515709 4.596724e+05  1735917.1
## 97   386.542191 4.694006e+05  1740016.1
## 98   390.568672 4.792308e+05  1752313.0
## 99   394.595153 4.891627e+05  1752313.0
## 100  398.621634 4.991966e+05  1776906.9
## 101  402.648115 5.093323e+05  1801500.8
## 102  406.674596 5.195699e+05  1887579.3
## 103  410.701078 5.299093e+05  1992103.3
## 104  414.727559 5.403506e+05  2004400.2
## 105  418.754040 5.508938e+05  2016697.1
## 106  422.780521 5.615389e+05  2005856.4
## 107  426.807002 5.722858e+05  2043232.6
## 108  430.833483 5.831345e+05  2049462.0
## 109  434.859964 5.940852e+05  2086838.2
## 110  438.886446 6.051377e+05  2111755.7
## 111  442.912927 6.162921e+05  2124214.4
## 112  446.939408 6.275483e+05  2149131.9
## 113  450.965889 6.389064e+05  2139912.4
## 114  454.992370 6.503664e+05  2168830.1
## 115  459.018851 6.619283e+05  2181625.6
## 116  463.045333 6.735920e+05  2220012.0
## 117  467.071814 6.853575e+05  2230967.0
## 118  471.098295 6.972250e+05  2250423.1
## 119  475.124776 7.091943e+05  2250423.1
## 120  479.151257 7.212655e+05  2276364.6
## 121  483.177738 7.334385e+05  2276364.6
## 122  487.204219 7.457134e+05  2347155.6
## 123  491.230701 7.580902e+05  2413836.1
## 124  495.257182 7.705688e+05  2453844.5
## 125  499.283663 7.831493e+05  2467180.6
## 126  503.310144 7.958317e+05  2441556.2
## 127  507.336625 8.086160e+05  2476941.0
## 128  511.363106 8.215021e+05  2518109.0
## 129  515.389588 8.344900e+05  2545554.4
## 130  519.416069 8.475799e+05  2586722.4
## 131  523.442550 8.607716e+05  2638686.9
## 132  527.469031 8.740652e+05  2663938.0
## 133  531.495512 8.874606e+05  2706334.9
## 134  535.521993 9.009579e+05  2734599.4
## 135  539.548474 9.145571e+05  2734599.4
## 136  543.574956 9.282581e+05  2734599.4
## 137  547.601437 9.420610e+05  2741665.6
## 138  551.627918 9.559658e+05  2826245.3
## 139  555.654399 9.699724e+05  2854938.1
## 140  559.680880 9.840809e+05  2912323.8
## 141  563.707361 9.982913e+05  2912323.8
## 142  567.733843 1.012604e+06  2941016.6
## 143  571.760324 1.027018e+06  2984055.9
## 144  575.786805 1.041534e+06  2984055.9
## 145  579.813286 1.056151e+06  2998402.3
## 146  583.839767 1.070871e+06  3012748.7
## 147  587.866248 1.085693e+06  3166076.1
## 148  591.892729 1.100616e+06  3173473.5
## 149  595.919211 1.115641e+06  3188268.3
## 150  599.945692 1.130769e+06  3210460.4
## 151  603.972173 1.145998e+06  3225255.1
## 152  607.998654 1.161329e+06  3225255.1
## 153  612.025135 1.176761e+06  3236411.8
## 154  616.051616 1.192296e+06  3251934.2
## 155  620.078097 1.207932e+06  3251934.2
## 156  624.104579 1.223671e+06  3267456.5
## 157  628.131060 1.239511e+06  3282978.9
## 158  632.157541 1.255453e+06  3298501.3
## 159  636.184022 1.271497e+06  3360590.7
## 160  640.210503 1.287643e+06  3376113.1
## 161  644.236984 1.303891e+06  3391635.4
## 162  648.263466 1.320240e+06  3391635.4
## 163  652.289947 1.336692e+06  3391635.4
## 164  656.316428 1.353245e+06  3407157.8
## 165  660.342909 1.369900e+06  3506606.1
## 166  664.369390 1.386657e+06  3538703.2
## 167  668.395871 1.403516e+06  3538703.2
## 168  672.422352 1.420477e+06  3570800.3
## 169  676.448834 1.437539e+06  3586848.8
## 170  680.475315 1.454704e+06  3618945.9
## 171  684.501796 1.471970e+06  3683140.1
## 172  688.528277 1.489339e+06  3683140.1
## 173  692.554758 1.506809e+06  3699188.6
## 174  696.581239 1.524381e+06  3731285.7
## 175  700.607721 1.542054e+06  3731285.7
## 176  704.634202 1.559830e+06  3753639.0
## 177  708.660683 1.577708e+06  3787455.6
## 178  712.687164 1.595687e+06  3804363.8
## 179  716.713645 1.613768e+06  3855088.7
## 180  720.740126 1.631952e+06  3863542.8
## 181  724.766607 1.650237e+06  3871997.0
## 182  728.793089 1.668624e+06  3905813.6
## 183  732.819570 1.687112e+06  3939630.1
## 184  736.846051 1.705703e+06  3964992.5
## 185  740.872532 1.724395e+06  3964992.5
## 186  744.899013 1.743190e+06  3976828.3
## 187  748.925494 1.762086e+06  3994044.0
## 188  752.951976 1.781084e+06  4080122.6
## 189  756.978457 1.800184e+06  4080122.6
## 190  761.004938 1.819386e+06  4080122.6
## 191  765.031419 1.838690e+06  4131769.7
## 192  769.057900 1.858095e+06  4209240.4
## 193  773.084381 1.877603e+06  4235063.9
## 194  777.110862 1.897212e+06  4235063.9
## 195  781.137344 1.916923e+06  4286711.1
## 196  785.163825 1.936736e+06  4286711.1
## 197  789.190306 1.956651e+06  4381397.5
## 198  793.216787 1.976668e+06  4424436.7
## 199  797.243268 1.996786e+06  4424436.7
## 200  801.269749 2.017007e+06  4510515.3
## 201  805.296231 2.037329e+06  4587985.9
## 202  809.322712 2.057753e+06  4639633.1
## 203  813.349193 2.078279e+06  4717103.7
## 204  817.375674 2.098907e+06  4734319.5
## 205  821.402155 2.119637e+06  4734319.5
## 206  825.428636 2.140469e+06  4734319.5
## 207  829.455117 2.161403e+06  4760621.2
## 208  833.481599 2.182438e+06  4760621.2
## 209  837.508080 2.203575e+06  4786923.0
## 210  841.534561 2.224814e+06  4821992.0
## 211  845.561042 2.246155e+06  4874595.6
## 212  849.587523 2.267598e+06  4909664.6
## 213  853.614004 2.289143e+06  4909664.6
## 214  857.640485 2.310790e+06  4975501.8
## 215  861.666967 2.332538e+06  4975501.8
## 216  865.693448 2.354389e+06  5029097.8
## 217  869.719929 2.376341e+06  5029097.8
## 218  873.746410 2.398395e+06  5046963.2
## 219  877.772891 2.420551e+06  5082693.9
## 220  881.799372 2.442809e+06  5134915.7
## 221  885.825854 2.465168e+06  5217034.4
## 222  889.852335 2.487630e+06  5263449.3
## 223  893.878816 2.510193e+06  5291298.2
## 224  897.905297 2.532859e+06  5328430.1
## 225  901.931778 2.555626e+06  5365562.0
## 226  905.958259 2.578495e+06  5384128.0
## 227  909.984740 2.601466e+06  5439825.9
## 228  914.011222 2.624538e+06  5523372.7
## 229  918.037703 2.647713e+06  5560504.6
## 230  922.064184 2.670989e+06  5634768.4
## 231  926.090665 2.694368e+06  5634768.4
## 232  930.117146 2.717848e+06  5681183.3
## 233  934.143627 2.741430e+06  5785338.4
## 234  938.170109 2.765114e+06  5823212.9
## 235  942.196590 2.788900e+06  5851618.8
## 236  946.223071 2.812788e+06  5889493.4
## 237  950.249552 2.836777e+06  5955773.9
## 238  954.276033 2.860869e+06  6096644.0
## 239  958.302514 2.885062e+06  6096644.0
## 240  962.328995 2.909357e+06  6096644.0
## 241  966.355477 2.933754e+06  6125629.7
## 242  970.381958 2.958253e+06  6135291.5
## 243  974.408439 2.982854e+06  6173939.0
## 244  978.434920 3.007556e+06  6318867.2
## 245  982.461401 3.032361e+06  6328529.1
## 246  986.487882 3.057267e+06  6328529.1
## 247  990.514364 3.082275e+06  6338190.9
## 248  994.540845 3.107385e+06  6357514.7
## 249  998.567326 3.132597e+06  6396162.2
## 250 1002.593807 3.157911e+06  6415486.0
## 251 1006.620288 3.183327e+06  6415486.0
## 252 1010.646769 3.208844e+06  6434809.7
## 253 1014.673250 3.234464e+06  6560414.1
## 254 1018.699732 3.260185e+06  6608723.5
## 255 1022.726213 3.286008e+06  6628047.2
## 256 1026.752694 3.311933e+06  6666694.7
## 257 1030.779175 3.337960e+06  6618184.1
## 258 1034.805656 3.364089e+06  6618184.1
## 259 1038.832137 3.390319e+06  6618184.1
## 260 1042.858619 3.416652e+06  6647773.6
## 261 1046.885100 3.443086e+06  6706952.6
## 262 1050.911581 3.469623e+06  6668339.3
## 263 1054.938062 3.496261e+06  6728777.4
## 264 1058.964543 3.523001e+06  6728777.4
## 265 1062.991024 3.549842e+06  6728777.4
## 266 1067.017505 3.576786e+06  6748923.5
## 267 1071.043987 3.603832e+06  6789215.6
## 268 1075.070468 3.630979e+06  6849653.7
## 269 1079.096949 3.658228e+06  6900018.8
## 270 1083.123430 3.685579e+06  7029435.2
## 271 1087.149911 3.713032e+06  7101479.2
## 272 1091.176392 3.740587e+06  7204399.2
## 273 1095.202873 3.768244e+06  7245567.2
## 274 1099.229355 3.796003e+06  7531871.9
## 275 1103.255836 3.823863e+06  7618951.3
## 276 1107.282317 3.851825e+06  7618951.3
## 277 1111.308798 3.879890e+06  7805853.5
## 278 1115.335279 3.908056e+06  7955966.1
## 279 1119.361760 3.936324e+06  7979060.4
## 280 1123.388242 3.964694e+06  8036796.0
## 281 1127.414723 3.993165e+06  8059890.2
## 282 1131.441204 4.021739e+06  8059890.2
## 283 1135.467685 4.050414e+06  8059890.2
## 284 1139.494166 4.079191e+06  8082984.4
## 285 1143.520647 4.108071e+06  8036203.8
## 286 1147.547128 4.137052e+06  8036203.8
## 287 1151.573610 4.166134e+06  8036203.8
## 288 1155.600091 4.195319e+06  8110636.8
## 289 1159.626572 4.224606e+06  8110636.8
## 290 1163.653053 4.253994e+06  8222765.4
## 291 1167.679534 4.283485e+06  8222765.4
## 292 1171.706015 4.313077e+06  8222765.4
## 293 1175.732497 4.342771e+06  8322435.3
## 294 1179.758978 4.372567e+06  8322435.3
## 295 1183.785459 4.402465e+06  8372270.2
## 296 1187.811940 4.432464e+06  8397187.7
## 297 1191.838421 4.462566e+06  8447022.6
## 298 1195.864902 4.492769e+06  8447022.6
## 299 1199.891383 4.523075e+06  8470593.2
## 300 1203.917865 4.553482e+06  8470593.2
## 301 1207.944346 4.583991e+06  8521775.0
## 302 1211.970827 4.614602e+06  8547365.9
## 303 1215.997308 4.645314e+06  8560161.4
## 304 1220.023789 4.676129e+06  8560161.4
## 305 1224.050270 4.707045e+06  8624138.7
## 306 1228.076752 4.738064e+06  8636934.1
## 307 1232.103233 4.769184e+06  8662525.1
## 308 1236.129714 4.800406e+06  8916301.6
## 309 1240.156195 4.831730e+06  8929452.5
## 310 1244.182676 4.863156e+06  8929452.5
## 311 1248.209157 4.894683e+06  8929452.5
## 312 1252.235638 4.926313e+06  8955754.3
## 313 1256.262120 4.958044e+06  8982056.1
## 314 1260.288601 4.989878e+06  9087263.2
## 315 1264.315082 5.021813e+06  9087263.2
## 316 1268.341563 5.053850e+06  9113565.0
## 317 1272.368044 5.085989e+06  9192470.3
## 318 1276.394525 5.118229e+06  9205621.2
## 319 1280.421006 5.150572e+06  9258224.7
## 320 1284.447488 5.183016e+06  9258224.7
## 321 1288.473969 5.215563e+06  9271375.6
## 322 1292.500450 5.248211e+06  9271375.6
## 323 1296.526931 5.280961e+06  9337130.0
## 324 1300.553412 5.313813e+06  9416035.4
## 325 1304.579893 5.346767e+06  9429186.2
## 326 1308.606375 5.379822e+06  9441585.7
## 327 1312.632856 5.412980e+06  9441585.7
## 328 1316.659337 5.446239e+06  9495692.2
## 329 1320.685818 5.479601e+06  9566110.2
## 330 1324.712299 5.513064e+06  9566110.2
## 331 1328.738780 5.546629e+06  9566110.2
## 332 1332.765261 5.580296e+06  9566110.2
## 333 1336.791743 5.614064e+06  9635732.5
## 334 1340.818224 5.647935e+06  9663581.5
## 335 1344.844705 5.681907e+06  9663581.5
## 336 1348.871186 5.715982e+06  9669488.8
## 337 1352.897667 5.750158e+06  9669488.8
## 338 1356.924148 5.784436e+06  9698181.7
## 339 1360.950630 5.818816e+06  9726874.5
## 340 1364.977111 5.853298e+06  9741220.9
## 341 1369.003592 5.887881e+06  9741220.9
## 342 1373.030073 5.922567e+06  9784260.2
## 343 1377.056554 5.957354e+06  9812953.1
## 344 1381.083035 5.992244e+06  9850438.9
## 345 1385.109516 6.027235e+06  9850438.9
## 346 1389.135998 6.062328e+06  9942070.9
## 347 1393.162479 6.097522e+06  9987886.9
## 348 1397.188960 6.132819e+06 10003158.8
## 349 1401.215441 6.168218e+06 10048974.8
## 350 1405.241922 6.203718e+06 10064246.8
## 351 1409.268403 6.239321e+06 10099881.5
## 352 1413.294885 6.275025e+06 10099881.5
## 353 1417.321366 6.310831e+06 10154299.0
## 354 1421.347847 6.346739e+06 10195695.1
## 355 1425.374328 6.382749e+06 10195695.1
## 356 1429.400809 6.418860e+06 10246420.0
## 357 1433.427290 6.455074e+06 10263328.2
## 358 1437.453771 6.491389e+06 10280236.5
## 359 1441.480253 6.527806e+06 10280236.5
## 360 1445.506734 6.564326e+06 10297144.8
## 361 1449.533215 6.600947e+06 10297144.8
## 362 1453.559696 6.637669e+06 10288040.4
## 363 1457.586177 6.674494e+06 10342667.1
## 364 1461.612658 6.711421e+06 10342667.1
## 365 1465.639140 6.748449e+06 10342667.1
## 366 1469.665621 6.785579e+06 10342667.1
## 367 1473.692102 6.822812e+06 10342667.1
## 368 1477.718583 6.860146e+06 10360876.0
## 369 1481.745064 6.897582e+06 10397293.9
## 370 1485.771545 6.935119e+06 10415502.8
## 371 1489.798026 6.972759e+06 10451920.6
## 372 1493.824508 7.010501e+06 10470129.6
## 373 1497.850989 7.048344e+06 10621342.8
## 374 1501.877470 7.086289e+06 10621342.8
## 375 1505.903951 7.124336e+06 10621342.8
## 376 1509.930432 7.162485e+06 10683094.8
## 377 1513.956913 7.200736e+06 10703678.8
## 378 1517.983394 7.239089e+06 10703678.8
## 379 1522.009876 7.277544e+06 10703678.8
## 380 1526.036357 7.316100e+06 10759816.9
## 381 1530.062838 7.354758e+06 10759816.9
## 382 1534.089319 7.393519e+06 10753668.5
## 383 1538.115800 7.432381e+06 10821301.6
## 384 1542.142281 7.471344e+06 10843846.0
## 385 1546.168763 7.510410e+06 10843846.0
## 386 1550.195244 7.549578e+06 10843846.0
## 387 1554.221725 7.588847e+06 10843846.0
## 388 1558.248206 7.628219e+06 10934023.5
## 389 1562.274687 7.667692e+06 10934023.5
## 390 1566.301168 7.707267e+06 10934023.5
## 391 1570.327649 7.746944e+06 10934023.5
## 392 1574.354131 7.786723e+06 11362366.7
## 393 1578.380612 7.826604e+06 11409709.9
## 394 1582.407093 7.866586e+06 11486954.0
## 395 1586.433574 7.906671e+06 11486954.0
## 396 1590.460055 7.946857e+06 11599082.7
## 397 1594.486536 7.987145e+06 11651686.2
## 398 1598.513018 8.027535e+06 11677988.0
## 399 1602.539499 8.068027e+06 11677988.0
## 400 1606.565980 8.108621e+06 11677988.0
## 401 1610.592461 8.149317e+06 11677988.0
## 402 1614.618942 8.190114e+06 11704289.8
## 403 1618.645423 8.231014e+06 11704289.8
## 404 1622.671904 8.272015e+06 11730591.5
## 405 1626.698386 8.313118e+06 11756893.3
## 406 1630.724867 8.354323e+06 11835798.6
## 407 1634.751348 8.395630e+06 11835798.6
## 408 1638.777829 8.437039e+06 11914704.0
## 409 1642.804310 8.478549e+06 11914704.0
## 410 1646.830791 8.520162e+06 11914704.0
## 411 1650.857273 8.561876e+06 11941005.7
## 412 1654.883754 8.603692e+06 11967307.5
## 413 1658.910235 8.645610e+06 11967307.5
## 414 1662.936716 8.687630e+06 11967307.5
## 415 1666.963197 8.729752e+06 11993609.3
## 416 1670.989678 8.771975e+06 11993609.3
## 417 1675.016159 8.814301e+06 11993609.3
## 418 1679.042641 8.856728e+06 11993609.3
## 419 1683.069122 8.899258e+06 11993609.3
## 420 1687.095603 8.941889e+06 11993609.3
## 421 1691.122084 8.984622e+06 12019911.1
## 422 1695.148565 9.027457e+06 12072514.6
## 423 1699.175046 9.070393e+06 12072514.6
## 424 1703.201528 9.113432e+06 12072514.6
## 425 1707.228009 9.156572e+06 12072514.6
## 426 1711.254490 9.199815e+06 12098816.4
## 427 1715.280971 9.243159e+06 12098816.4
## 428 1719.307452 9.286605e+06 12125118.2
## 429 1723.333933 9.330153e+06 12125118.2
## 430 1727.360414 9.373802e+06 12125118.2
## 431 1731.386896 9.417554e+06 12125118.2
## 432 1735.413377 9.461408e+06 12125118.2
## 433 1739.439858 9.505363e+06 12125118.2
## 434 1743.466339 9.549420e+06 12125118.2
## 435 1747.492820 9.593579e+06 12125118.2
## 436 1751.519301 9.637840e+06 12125118.2
## 437 1755.545782 9.682203e+06 12125118.2
## 438 1759.572264 9.726668e+06 12125118.2
## 439 1763.598745 9.771234e+06 12125118.2
## 440 1767.625226 9.815903e+06 12151419.9
## 441 1771.651707 9.860673e+06 12151419.9
## 442 1775.678188 9.905545e+06 12177721.7
## 443 1779.704669 9.950519e+06 12177721.7
## 444 1783.731151 9.995595e+06 12058590.1
## 445 1787.757632 1.004077e+07 12058590.1
## 446 1791.784113 1.008605e+07 12058590.1
## 447 1795.810594 1.013143e+07 12058590.1
## 448 1799.837075 1.017692e+07 12058590.1
## 449 1803.863556 1.022250e+07 12058590.1
## 450 1807.890037 1.026819e+07 12058590.1
## 451 1811.916519 1.031398e+07 12058590.1
## 452 1815.943000 1.035987e+07 12058590.1
## 453 1819.969481 1.040586e+07 12058590.1
## 454 1823.995962 1.045196e+07 12114288.0
## 455 1828.022443 1.049815e+07 12197834.8
## 456 1832.048924 1.054445e+07 12197834.8
## 457 1836.075406 1.059085e+07 12197834.8
## 458 1840.101887 1.063735e+07 12197834.8
## 459 1844.128368 1.068396e+07 12197834.8
## 460 1848.154849 1.073066e+07 12225683.8
## 461 1852.181330 1.077747e+07 12225683.8
## 462 1856.207811 1.082438e+07 12225683.8
## 463 1860.234292 1.087139e+07 12225683.8
## 464 1864.260774 1.091851e+07 12225683.8
## 465 1868.287255 1.096572e+07 12250051.6
## 466 1872.313736 1.101304e+07 12250051.6
## 467 1876.340217 1.106046e+07 12279641.1
## 468 1880.366698 1.110798e+07 12279641.1
## 469 1884.393179 1.115560e+07 12279641.1
## 470 1888.419661 1.120332e+07 12338820.1
## 471 1892.446142 1.125115e+07 12467041.2
## 472 1896.472623 1.129908e+07 12530165.5
## 473 1900.499104 1.134711e+07 12530165.5
## 474 1904.525585 1.139524e+07 12530165.5
## 475 1908.552066 1.144347e+07 12512130.0
## 476 1912.578547 1.149181e+07 12545946.6
## 477 1916.605029 1.154025e+07 12579763.1
## 478 1920.631510 1.158879e+07 12613579.7
## 479 1924.657991 1.163743e+07 12613579.7
## 480 1928.684472 1.168617e+07 12681212.8
## 481 1932.710953 1.173502e+07 12681212.8
## 482 1936.737434 1.178396e+07 12681212.8
## 483 1940.763916 1.183301e+07 12715029.4
## 484 1944.790397 1.188216e+07 12673409.0
## 485 1948.816878 1.193141e+07 12709826.8
## 486 1952.843359 1.198077e+07 12709826.8
## 487 1956.869840 1.203023e+07 12709826.8
## 488 1960.896321 1.207978e+07 12782662.5
## 489 1964.922802 1.212944e+07 12610505.5
## 490 1968.949284 1.217920e+07 12610505.5
## 491 1972.975765 1.222907e+07 12610505.5
## 492 1977.002246 1.227903e+07 12545946.6
## 493 1981.028727 1.232910e+07 12545946.6
## 494 1985.055208 1.237927e+07 12545946.6
## 495 1989.081689 1.242954e+07 12593289.8
## 496 1993.108170 1.247991e+07 12593289.8
## 497 1997.134652 1.253039e+07 12593289.8
## 498 2001.161133 1.258097e+07 12593289.8
## 499 2005.187614 1.263164e+07 12677455.4
## 500 2009.214095 1.268243e+07 12677455.4
## 501 2013.240576 1.273331e+07 12677455.4
## 502 2017.267057 1.278429e+07 12677455.4
## 503 2021.293539 1.283538e+07 12677455.4
## 504 2025.320020 1.288657e+07 12677455.4
## 505 2029.346501 1.293786e+07 12677455.4
## 506 2033.372982 1.298925e+07 12677455.4
## 507 2037.399463 1.304074e+07 12677455.4
## 508 2041.425944 1.309234e+07 12677455.4
## 509 2045.452425 1.314403e+07 12730059.0
## 510 2049.478907 1.319583e+07 12730059.0
## 511 2053.505388 1.324773e+07 12730059.0
## 512 2057.531869 1.329974e+07 12730059.0
## 513 2061.558350 1.335184e+07 12730059.0

There are a lot of elements i the plot of K that can be explained. First, the Kpois(r) line in Red is the theoretical value of K for each distance window (r) under a Poisson assumption of Complete Spatial Randomness.[Practical book] The Black line is the estimated values of K accounting for the effects of the edge of the study area.[Practical book]

(Here, the correction specifies how points towards the edge are dealt with, in this case, border means that points towards the edge are ignored for the calculation but are included for the central points.)

If the value of K falls above the line, it means that there is cluster of the data at the distance. If the value of K is below the line, there shows a dispersion of the data. From the graph, we can see that up until distances of around 2000 metres, the location of the eviction seems to be clustered in Harrow. (however, at around 1500 m, the distribution appears random and then dispersed between about 1600 and 2100 metres.)

Density-based spatial clustering of applications with noise: DBSCAN

# Load the package
library(raster)

## 
## Attaching package: 'raster'

## The following object is masked from 'package:janitor':
## 
##     crosstab

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:tidyr':
## 
##     extract

## The following objects are masked from 'package:MASS':
## 
##     area, select

## The following object is masked from 'package:nlme':
## 
##     getData

## The following objects are masked from 'package:spatstat.geom':
## 
##     area, rotate, shift

library(fpc)

#first check the coordinate reference system of the CD502 spatial polygon:
st_geometry(NYDistricts)

## Geometry set for 71 features 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 563069.7 ymin: 4483098 xmax: 609762.3 ymax: 4529952
## Projected CRS: WGS 84 / UTM zone 18N
## First 5 geometries:

## MULTIPOLYGON (((595104.1 4522026, 595098.5 4521...

## MULTIPOLYGON (((597308.6 4510424, 597266.4 4510...

## MULTIPOLYGON (((594364.8 4521323, 594363 452132...

## MULTIPOLYGON (((593253.9 4504080, 593249.3 4504...

## MULTIPOLYGON (((593432.1 4523996, 593429 452399...

DBSCAN requires to input two parameters: 1. Epsilon - this is the radius within which the algorithm with search for clusters 2. MinPts - this is the minimum number of points that should be considered a cluster

Based on the results of the Ripley’s K analysis earlier, we can see that we are getting clustering up to a radius of around 2000m, with the largest bulge in the graph at around 1700m. Therefore, 1700m is probably a good place to start and we will begin by searching for clusters of at least 4 points…

#first extract the points from the spatial points data frame
NewEvictionsCDPoints <- NewEvictionsCD %>%
  coordinates(.)%>%
  as.data.frame()

#now run the dbscan analysis
db <- NewEvictionsCDPoints %>%
  fpc::dbscan(.,eps = 800, MinPts = 6)

#now plot the result
plot(db, NewEvictionsCDPoints, main = "DBSCAN Output", frame = F)
plot(NYDistricts$geometry, add=T)

Use kNNdistplot() to find a suitable eps value based on the ‘knee’ in the plot…

# k is no of nearest neighbours used, use min points
library(dbscan)

## 
## Attaching package: 'dbscan'

## The following object is masked from 'package:fpc':
## 
##     dbscan

NewEvictionsCDPoints%>%
  dbscan::kNNdistplot(.,k=6)

https://www.datanovia.com/en/lessons/dbscan-density-based-clustering-essentials/

Advanced DBSCAN

library(ggplot2)

db

## dbscan Pts=116 MinPts=6 eps=800
##         0  1  2  3 4
## border 22  3  1  4 0
## seed    0 21 25 32 8
## total  22 24 26 36 8

db$cluster

##   [1] 0 3 0 1 1 0 2 3 0 3 2 2 3 2 3 4 3 3 3 3 1 3 3 3 2 2 3 3 0 3 1 1 4 1 0 3 3
##  [38] 0 2 2 2 0 2 3 2 0 1 0 3 3 3 1 3 2 3 3 2 3 3 0 2 1 1 1 0 3 2 2 2 4 3 1 1 0
##  [75] 2 4 0 1 2 0 0 3 1 2 3 1 4 3 4 2 1 1 3 1 0 3 2 0 1 2 2 1 0 4 1 0 0 3 3 4 2
## [112] 3 1 2 0 1

# add this cluster membership info back into our dataframe
NewEvictionsCDPoints <- NewEvictionsCDPoints %>%
  mutate(dbcluster=db$cluster)

# Next we are going to create some convex hull polygons to wrap around the points in our clusters.
chulls <- NewEvictionsCDPoints %>%
  group_by(dbcluster) %>%
  dplyr::mutate(hull = 1:n(),
  hull = factor(hull, chull(coords.x1, coords.x2)))%>%
  arrange(hull)

#chulls2 <- ddply(BluePlaquesSubPoints, .(dbcluster), 
              #  function(df) df[chull(df$coords.x1, df$coords.x2), ])

# As 0 isn’t actually a cluster (it’s all points that aren’t in a cluster) drop it from the dataframe
chulls <- chulls %>%
  filter(dbcluster >=1)

# ggplot
dbplot <- ggplot(data = NewEvictionsCDPoints, 
                 aes(coords.x1,coords.x2, colour=dbcluster, fill=dbcluster)) 
#add the points in
dbplot <- dbplot + geom_point()
#now the convex hulls
dbplot <- dbplot + geom_polygon(data = chulls, 
                                aes(coords.x1,coords.x2, group=dbcluster), 
                                alpha = 0.5) 
#now plot, setting the coordinates to scale correctly and as a black and white plot 
#(just for the hell of it)...
dbplot + theme_bw() + coord_equal()

Nicer way of DBSCAN

#convex hulls to wrap around points
chulls2 <- data.frame()
for (cluster in 1:max(NewEvictionsCDPoints$dbcluster)) {
  cluster_data <- NewEvictionsCDPoints %>%
    filter(dbcluster == cluster)
  ch <- chull(cluster_data$coords.x1, cluster_data$coords.x2)
  chulls2 <- chulls2 %>%
    bind_rows(cluster_data[c(ch), ])
}

# ggplot
dbplot <- ggplot(data = NewEvictionsCDPoints, 
                 aes(coords.x1,coords.x2, colour=dbcluster, fill=dbcluster)) 
#add the points in
dbplot <- dbplot + geom_point()
#now the convex hulls
dbplot <- dbplot + geom_polygon(data = chulls2, 
                                aes(coords.x1,coords.x2, group=dbcluster), 
                                alpha = 0.5) 
#now plot, setting the coordinates to scale correctly and as a black and white plot 
#(just for the hell of it)...
dbplot + theme_bw() + coord_equal()

Moran’s I and Spatial Autocorrelation

library(here)
library(janitor)
library(dplyr)

# join the New York District data and eviction data
#Joinfun <- function(data1, data2){

#output<- data1%>%
#  st_join(NYDistricts,.)%>%
#   add_count(boro_cd, name="eviction-borough") 

#  return(output)
#}

#NYEvictions <- Joinfun(NewEvictions, NYDistricts)

#NYEvictions

#summary(NYEvictions)

#NewEvictions <- NewEvictions[NYEvictions,]

#tm_shape(NYDistricts) +
#  tm_polygons(col = NA, alpha = 0.5) +
#tm_shape(NYEvictions) +
#  tm_dots(col = "blue")

library(sf)
points_sf_joined <- NYDistricts %>%
  st_join(NewEvictions) %>%
  add_count(boro_cd) %>%
  janitor::clean_names() %>%
  #calculate area
  mutate(area=st_area(.)) %>%
  #then density of the points per ward
  mutate(density=n/area) %>%
  #select density and some other variables 
  dplyr::select(density, boro_cd, n)

points_sf_joined <- points_sf_joined %>%                    
  group_by(boro_cd) %>%         
  summarise(density = first(density),
          boro_cd = first(boro_cd),
          plaquecount= first(n))

tm_shape(points_sf_joined) +
    tm_polygons("density",
        style="jenks",
        palette="PuOr",
        midpoint=NA,
        popup.vars=c("boro_cd", "density"),
        title="Evictions Density")

So, from the map, it looks as though we might have some clustering of blue plaques in the Manhatten so let’s check this with Moran’s I and some other statistics.

Before being able to calculate Moran’s I and any similar statistics, we need to first define a W_ij spatial weights matrix.

library(spdep)

## Loading required package: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

#First calculate the centroids of all Wards in New York City

coordsW <- points_sf_joined%>%
  st_centroid()%>%
  st_geometry()

## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x

plot(coordsW,axes=TRUE)

#create a neighbours list
LWard_nb <- points_sf_joined %>%
  poly2nb(., queen=T)

summary(LWard_nb)

## Neighbour list object:
## Number of regions: 71 
## Number of nonzero links: 314 
## Percentage nonzero weights: 6.228923 
## Average number of links: 4.422535 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8 
##  1  8 11 16 18 10  6  1 
## 1 least connected region:
## 62 with 1 link
## 1 most connected region:
## 33 with 8 links

Average number of links is 4.422535

#plot the neighbours
plot(LWard_nb, st_geometry(coordsW), col="red")
#add a map underneath
plot(points_sf_joined$geometry, add=T)

#create a spatial weights matrix from these weights
#Lward.lw <- LWard_nb %>%
#  nb2mat(., style="B")

#sum(Lward.lw)
#summary(Lward.lw)

Spatial autocorrelation

Moran’s I

Now we have defined our W_ij matrix, we can calculate the Moran’s I and other associated statistics. However, Moran’s I requires a spatial weight list type object as opposed to matrix.

Lward.lw <- LWard_nb %>%
  nb2listw(., style="C")

Moran’s I test tells us whether we have clustered values (close to 1) or dispersed values (close to -1)

I_LWard_Global_Density <- points_sf_joined %>%
  pull(density) %>%
  as.vector()%>%
  moran.test(., Lward.lw)

I_LWard_Global_Density

## 
##  Moran I test under randomisation
## 
## data:  .  
## weights: Lward.lw    
## 
## Moran I statistic standard deviate = 8.238, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.593759367      -0.014285714       0.005447915

Geary’s C

This tells us whether similar values or dissimilar values are clustering

C_LWard_Global_Density <- 
  points_sf_joined %>%
  pull(density) %>%
  as.vector()%>%
  geary.test(., Lward.lw)

C_LWard_Global_Density

## 
##  Geary C test under randomisation
## 
## data:  . 
## weights: Lward.lw 
## 
## Geary C statistic standard deviate = 4.0732, p-value = 2.319e-05
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.46436106        1.00000000        0.01729311

Getis Ord General G

This tells us whether high or low values are clustering. If G > Expected = High values clustering; if G < expected = low values clustering

G_LWard_Global_Density <- 
  points_sf_joined %>%
  pull(density) %>%
  as.vector()%>%
  globalG.test(., Lward.lw)

G_LWard_Global_Density

## 
##  Getis-Ord global G statistic
## 
## data:  . 
## weights: Lward.lw 
## 
## standard deviate = 7.2857, p-value = 1.6e-13
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       3.048722e-02       1.428571e-02       4.945053e-06

So the global statistics are indicating that we have spatial autocorrelation of Blue Plaques in London:

The Moran’s I statistic = 0.59 (remember 1 = clustered, 0 = no pattern, -1 = dispersed) which shows that we have some distinctive clustering

The Geary’s C statistic = 0.46 (remember Geary’s C falls between 0 and 2; 1 means no spatial autocorrelation, <1 - positive spatial autocorrelation or similar values clustering, >1 - negative spatial autocorreation or dissimilar values clustering) which shows that similar values are clustering

The General G statistic = G > expected, so high values are tending to cluster.

We can now also calculate local versions of the Moran’s I statistic (for each Ward) and a Getis Ord G*_i statistic to see where we have hot-spots…

#use the localmoran function to generate I for each ward in the city

I_LWard_Local_count <- points_sf_joined %>%
  pull(plaquecount) %>%
  as.vector()%>%
  localmoran(., Lward.lw)%>%
  as_tibble()

I_LWard_Local_Density <- points_sf_joined %>%
  pull(density) %>%
  as.vector()%>%
  localmoran(., Lward.lw)%>%
  as_tibble()

#what does the output (the localMoran object) look like?
slice_head(I_LWard_Local_Density, n=5)

## # A tibble: 5 × 5
##   Ii           E.Ii          Var.Ii     Z.Ii         `Pr(z != E(Ii))`
##   <localmrn>   <localmrn>    <localmrn> <localmrn>   <localmrn>      
## 1 -0.001116426 -0.0015044943 0.02372412  0.002519496 0.9979897       
## 2 -0.007248032 -0.0007965765 0.01222156 -0.058357190 0.9534641       
## 3 -0.056804476 -0.0008635776 0.01324855 -0.486010031 0.6269601       
## 4 -0.158916086 -0.0041275381 0.06309378 -0.616233868 0.5377402       
## 5  0.011812369 -0.0008714819 0.01296872  0.111378795 0.9113160

points_sf_joined <- points_sf_joined %>%
  mutate(plaque_count_I = as.numeric(I_LWard_Local_count$Ii))%>%
  mutate(plaque_count_Iz =as.numeric(I_LWard_Local_count$Z.Ii))%>%
  mutate(density_I =as.numeric(I_LWard_Local_Density$Ii))%>%
  mutate(density_Iz =as.numeric(I_LWard_Local_Density$Z.Ii))

Mapping outputs

We’ll set the breaks manually based on the rule that data points >2.58 or <-2.58 standard deviations away from the mean are significant at the 99% level (<1% chance that autocorrelation not present); >1.96 - <2.58 or <-1.96 to >-2.58 standard deviations are significant at the 95% level (<5% change that autocorrelation not present).

breaks1<-c(-1000,-2.58,-1.96,-1.65,1.65,1.96,2.58,1000)

MoranColours<- rev(brewer.pal(8, "RdGy"))

tm_shape(points_sf_joined) +
    tm_polygons("plaque_count_Iz",
        style="fixed",
        breaks=breaks1,
        palette=MoranColours,
        midpoint=NA,
        title="Local Moran's I, Blue Plaques in London")

This map shows some areas in the centre of London that have relatively high scores, indicating areas with lots of blue plaques neighbouring other areas with lots of blue plaques.

the Getis Ord G*_i statisic for hot and cold spots

Gi_LWard_Local_Density <- points_sf_joined %>%
  pull(density) %>%
  as.vector()%>%
  localG(., Lward.lw)

head(Gi_LWard_Local_Density)

## [1] -0.002519496  0.058357190 -0.486010031 -0.616233868 -0.111378795
## [6] -0.479051443

points_sf_joined <- points_sf_joined %>%
  mutate(density_G = as.numeric(Gi_LWard_Local_Density))

GIColours<- rev(brewer.pal(8, "RdBu"))

#now plot on an interactive map
tm_shape(points_sf_joined) +
    tm_polygons("density_G",
        style="fixed",
        breaks=breaks1,
        palette=GIColours,
        midpoint=NA,
        title="Gi*, Blue Plaques in London")

#### Other variables can be checked in the csv file.

Other Analysis Methods